/*
Este programa implementa o algoritmo de Monte Carlo, assim como montecarlo.c.
Este, porém, utiliza Pthreads para paralelizar parte do código;
especificamente, cria tantas threads quanto especificado em NUM_THREADS, e roda
o loop principal em paralelo (como são amostragens independentes, é possível
paralelizar mais do que nos outros algoritmos).
*/

#include <pthread.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <gmp.h>
#define NUM_THREADS 4

//Valores totais de m e n (soma de todas as threads)
mpz_t total_m, total_n;

//Número de iterações totais e por thread (caso a divisão não seja exata, a
//última thread roda thread_remainder iterações
mpz_t total_iterations, thread_iterations, thread_remainder;

//Mutex para escrita a total_m e total_n
pthread_mutex_t mutex;

void* montecarlo(void *thread_id)
{
   double x, y;
   
   //A id da thread servirá como seed para os números pseudo-aleatórios
   int id;
   mpz_t m, n, *iterations;
   
   id = *((int *) thread_id);
   
   //Inicialização das variáveis
   mpz_init(m);
   mpz_init(n);
   
   //Verifica se é a última thread. Caso positivo, o número de iterações é
   //thread_remainder
   if(id == NUM_THREADS - 1)
      iterations = &thread_remainder;
   else
      iterations = &thread_iterations;

   //Loop principal
   while(mpz_cmp(n, *iterations) != 0)
   {
      //Amostra um ponto no quadrado unitário
      x = ((double) rand_r(&id) / (double) RAND_MAX);
      y = ((double) rand_r(&id) / (double) RAND_MAX);
      
      //Incrementa o número de iterações
      mpz_add_ui(n, n, 1);

      //Se amostra estiver dentro do círculo unitário, incrementa m
      if(sqrt(x*x + y*y) < 1.0)
         mpz_add_ui(m, m, 1);
   }   
   
   //Soma m e n aos m e n totais
   pthread_mutex_lock(&mutex);
   mpz_add(total_m, total_m, m);
   mpz_add(total_n, total_n, n);
   pthread_mutex_unlock(&mutex);
   
   pthread_exit(NULL);
}

int main(int argc, const char* argv[])
{
   //Declaração de variáveis para abertura de arquivo e medida de tempo
   FILE *fp;
   char filename[50];
   struct timeval ti, tf;
   long int sec, usec;

   //Início da contagem de tempo
   gettimeofday(&ti, NULL);

   //Declaração das threads e variáveis auxiliares
   pthread_t threads[NUM_THREADS];
   mpf_t pi, aux;
   int precision, t, rc;
   int thread_id[NUM_THREADS];
   void *status;

   //Inicialização de variáveis com a precisão especificada
   precision = atoi(argv[1]);
   mpf_set_default_prec(precision);
   mpz_init(total_m);
   mpz_init(total_n);
   mpz_init(total_iterations);
   mpz_init(thread_iterations);
   mpz_init(thread_remainder);
   mpf_init(pi);
   mpf_init(aux);
   mpz_ui_pow_ui(total_iterations, 10, precision + 1);

   //Inicializa mutex   
   pthread_mutex_init(&mutex, NULL);


   //Divisão das iterações entre as threads 
   mpz_cdiv_qr_ui(thread_iterations, thread_remainder, total_iterations, NUM_THREADS);
   if(mpz_cmp_ui(thread_remainder, 0) != 0)
      mpz_cdiv_qr_ui(thread_iterations, thread_remainder, total_iterations, NUM_THREADS - 1);
   else
      mpz_set(thread_remainder, thread_iterations);
      
   //Cria threads
   for(t = 0; t < NUM_THREADS; t++)
   {
      thread_id[t] = t;
      rc = pthread_create(&threads[t], NULL, montecarlo, (void *) &thread_id[t]);
      if (rc)
      {
         printf("ERRO; codigo  da  pthread_create() e %d\n", rc);
         exit(-1);
      }
   }
   
   //Aguarda as threads finalizarem
   for(t = 0; t < NUM_THREADS; t++)
      pthread_join(threads[t], &status);

   //Cálculo final de pi (pi =~ 4*m/n)
   mpf_set_z(pi, total_m);
   mpf_mul_ui(pi, pi, 4);
   mpf_set_z(aux, total_n);
   mpf_div(pi, pi, aux);
   
   
   //Término da contagem de tempo
   gettimeofday(&tf, NULL);
   sec = tf.tv_sec - ti.tv_sec;
   if(tf.tv_usec < ti.tv_usec)
      sec--;
   usec = (tf.tv_usec - ti.tv_usec + 1000000) % 1000000;
   printf("Tempo decorrido: %ld.%06ld s. Iterações: ", sec, usec);
   gmp_printf("%Zd\n", total_iterations);
   
   //Salva o resultado para arquivo
   sprintf(filename, "pi/montecarlo-seq-%d", precision);
   fp = fopen(filename, "w");
   gmp_fprintf(fp, "%.Ff\n", pi);
   fclose(fp);
   
   return 0;
}
